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ABSTRACT 

Imaging observations are generally affected by a fluctuating background of speckles, a particular problem 
when detecting faint stellar companions at small angular separations. These speckles can be created by both 
short-lived atmopsheric aberrations and slowly changing distortions in the optical system. Over the course of 
a long-exposure image, the combination of many independent realizations of speckle patterns forms a halo 
in the point spread function (PSF) of characteristic scale ^ A/ro (where ro is the coherence length in the 
pupil). While adaptive optics can increase the achievable image contrast, speckle noise remains a major source 
of random error, which decreases the sensitivity of companion detection observations near the diffraction limit. 
Knowing the distribution of the speckle intensities at a given location in the image plane is therefore important 
for understanding the noise limits of companion detection. 

The speckle noise limit in a long-exposure image is characterized by the intensity variance and the speckle 
lifetime. In this paper we address the former quantity through the distribution function of speckle intensity. 
Previous theoretical work has predicted a form for this distribution function at a single location in the image 
plane. We developed a fast readout mode to take short exposures of stellar images corrected by adaptive optics 
at the ground-based UCO/Lick Observatory, with integration times of 5 ms and a time between successive 
frames of 14.5 ms (A = 2.2 fim). These observations temporally oversample and spatially Nyquist sample the 
observed speckle patterns. We show, for various locations in the image plane, the observed distribution of 
speckle intensities is consistent with the predicted form. Additionally, we demonstrate a method by which Ic 
and Is can be mapped over the image plane. As the quantity Ic is proportional to the PSF of the telescope free 
of random atmospheric aberrations, this method can be used for PSF calibration and reconstruction. 
Subject headings: atmospheric effects — instrumentation: adaptive optics — techniques: high angular resolu- 
tion — techniques: image processing 



1. INTRODUCTION 

Wavefronts from astronomical sources are invariably dis- 
torted before being imaged by the observer, resulting spatio- 
temporally fluctuating "speckles" of image plane intensity. 
Speckle noise limits the dynami c range of long exposures in 
the v icinity of a bright source jRoddiedll981t iRacine et alJ 
11999*). hampering observational goals like direct detection of 
close-in extrasolar planetar y systems . 

Speckle interferometry ( lLabevriji 19701) and adaptive op- 
tics (AO) have succeeded in achieving diffraction-limited 
spatial resolution in the presence of atmospheric turbu- 
lence. However, contrast levels achievable with these 
techniques are limited by speckles and are generally or- 
ders of magnitude below the dynamic range required for 
planet detection. Additional strategies exist to remove 
speckle noise and increase contrast. If the speckles are 
long-lived as in space telescopes, reference PSF subtrac- 
tion and roll subtraction can be used (Fraauelli et al. 2004). 
More advanced techniques exist for speckles varying on 
shorter timescales, such as multiwavelength Simultaneous 
Differential I maging (Racine et al. 1999; Marois et al. 2000, 
I2004t 12005* 'Lenzen et al. 2004; Bifler et al. 2OO40, "dark" 
speckle (Labevrie 1995; Boccaletti et al. 1998), and Syn - 
chronous Interferometric Speckle Subtraction (Guvon 2004). 

We study the fundamental properties of speckles to illu- 
minate both the limits and design of future instrumentation 
and observational methodology, in particular high-order adap- 
tive optics coronagraphy (e .g. Sivaramakrishnan et al. 20Qil 
iFord et alJ20"03tlTrauger et all2003HMacintosh et alJ2004r 

' fitz@astron.berkeley.edu 



Several authors have developed statistical descriptions of 
stellar speckles in ground-based images. In the uncompen- 
sated case of the visible and near-IR, speckle interferometry is 
the main driver of these theoretical and observational studies. 
This method uses sequences of short-exposure, spatially well- 
sampled images to retrieve a single high-resolution image. 
The focus of previous statistical studies lies in short-exposure 
images' spatial power spectra and atmosphere/telescope opti- 
cal transfer functions which are central to the technique (e.g. 
Fried 1966; Daintv 1974; Daintv et al. 1981; Roddiet 1981] 
lATmee t al. 1986; Vernin et al. 1991). 

In contrast to techniques which increase resolution of un- 
compensated images through post-processing, AO images 
possess a diffraction-limited core formed by spatially coher- 
ent light in the pupil. The superposition of this coherent 
light with residual wave aberrations modifies the character of 
the speckle distribution from the uncompensated case. The 
description of speckle statistics in bot h regimes can be re - 
lated to the study of laser speckles (Goodmanl l 197511 198.^ . 
the theoretical statistical formalisms of which have been pre- 
viousl y applied to AO-corrected astron omical sp eckle pat- 
terns (Canales & Cagigal 1222^; Soummer & Aims .20041 
l^me & Soummer 2004, hereafter AS04). 

For companion detection, we are interested in the dis- 
tribution of intensity in the image plane. In high-Strehl- 
ratio images, diffraction of wave aberrations results in 
amplification of speckles near Airy maxima, hence the 
phenomenon of "speckle pinning" ( Bloemhof et al. 200 jl 
i Sivaramakrishnan et a l.' '2002', 'Bloemhof i2003i: .Perrin et all 
I2003h . As noted by AS04, this behavior is described 
by the theoretical distribution of intensity originally de- 
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rived by iGoodmanI (1 19751) (see also iGoodmanl 119851) . 
ICagigal & Canalesl (12001) observed this distribution in a lab- 
oratory setting, using a spatial light modulator to simultane- 
ously simulate turbulence and partial correction with AO. 

This paper presents on-sky results consistent with predic- 
tions of the theoretical speckle statistics for compensated as- 
tronomical imaging through the atmosphere. Section Impar- 
tially recapitulates the treatment of the intensity distribution 
given by AS04. Section |3 describes observations of short- 
exposure stellar images under moderate correction. Analysis 
and discussion of these results in light of direct imaging tech- 
niques follow in sections|4]and|5] 

2. PROBABILISTIC DESCRIPTION OF INTENSITY 

Here we characterize the probability distributions of im- 
age plane intensity in the case of an unresolved astronomi- 
cal source subject to random wavefront aberrations. In this 
section, we follow the development of the probability den- 
sity function presented by AS04, adopting their notation. 
The problem has been studied elsewhere to the same conclu- 
sions^by Goodman ( 1975, 1985). Canales & CagigaU1999b.) . 
and iSoummer & Aimd (I2004f) . 

2.1. Intensity in the Image Plane 

We wish to study the properties of intensity in the image 
plane by considering an undistorted light wave from an astro- 
nomical source, as well as aberrations of the wave introduced 
by the atmosphere and telescope. Studying the superposition 
of the perfect wave and the aberrations allows us to conve- 
niently express intensity statistics useful for companion de- 
tection. Here, we consider the case of a narrow-band wave 
(AA/A <C 1) which has entered the telescope, passed through 
any corrective optics (e.g. AO system), and further propa- 
gated into the entrance pupil of an imaging system. We ignore 
polarization and treat the light as a scalar wave. We represent 
the wave at a given location by a complex amplitude, which 
describes the envelope function of the oscillating (scalar) elec- 
tric field. 

At a given time f and position u in the entrance pupil of the 
imaging system, the complex amplitude of the wave is repre- 
sented by ^'i(u,f) and can be decomposed into the superposi- 
tion of an undisturbed wave of amplitude A, which is constant 
in space and time, and random wave aberrations a(u, t), 

*,(u,f)=[A + fl(u,f)]nu). (1) 

The aberrations a{u,t) are time-varying fluctuations in the 
residual wave amplitude arising from the partially-corrected 
distortions. The pupil function (unity inside the pupil, zero 
outside) and static aberrations are encapsulated in P(u). With 
this definition we may take A as a positive real number with- 
out loss of generality. 

The wave in the image plane ^'2(x) is the Fourier transform 
^{•} of the wave in the pupil plane, 

^'2(x, f ) = * 1 (u, f)} = AJ^{P(u)} + J^{a(u, t)P(u)}, 

= Cix) + S{x,t). (2) 

The resulting image plane wave is represented as a superposi- 
tion of two complex wave amplitudes. C represents the ampli- 
tude of the temporally static portion of the image plane wave, 
while S represents the amplitude arising from random aberra- 
tions. 

The intensity is the squared modulus of ^'2(x,f), 

/(x,f)= |*2(x,f)|'= |C(x)|2+|5(x,f)|' + 2Re[C*(x)5(x,f)] , 

(3) 



where * represents the complex conjugate. Additionally, one 
can define and /, as the time-averaged intensities corre- 
sponding to the coherent and speckle wavefront amplitudes 
in the image plane: 

(/(x,f))=/c(x)+/.(x) (4) 
/,(x)=|C(x)|2 (5) 

/,(x)=(|5(x,f)|') (6) 
The final intensity term on the right of equation Q is an in- 
terference term resulting from the superposition of the C and 
S waves, which has zero mean, and thus does not enter into 
equation 0. 

Although Ic{x) and /j(x) are positive definitequantities, we 
stress that energy is conserved when equation (|4j is integrated 
over the image plane. Indeed, when adding perturbations to 
the system, energy conservation constrains the value of A. 
Consider the case of pure phase perturbations, 

*i(u,f) = e''^*"-'*i'(u), (7) 
which have no effect on the spatially-integrated mean inten- 
sity. Qualitatively, as a direct result of the definition in equa- 
tion Q and energy conservation, an increase in the scale of 
0(u,f) decreases the value of A and the scale of J Idxjdx. 

In this paper, probabilistic descriptions are formed of the in- 
tensity at a single location in space, /(?), omitting x for brevity. 
We describe the intensity as a stochastic process - a sample 
function /(f) is observed, which is drawn from the overall en- 
semble of such functions. For simplicity, we assume that the 
aberrations driving the intensity fluctuations constitute an er- 
godic process, so that we may take ensemble averages (•) in 
equations (0} and (|6} to equivalently represent time averages. 
In reality, atmospheric fluctuations are not temporally station- 
ary, and we could redefine A as a function of time, and Ic, Is 
as the ensemble-averages of temporal functions. 

2.2. The Distribution of Intensity 

The behavior of speckles can be described via the probabil- 
ity distribution of intensity, characterized by /?(/), the proba- 
bility density function {pdf). This is a first-order description 
of the stochastic process for / in which time is marginalized. 

We assume that in equation Q, the image-plane amplitude 
'J'2 is formed by the linear combination of a large number of 
independent phasors from the pupil plane. From the central 
limit theorem, this addition results in Gaussian statistics for 
the real and imaginary parts of the image plane amplitude, 
regardless of the particular statistics governing the pupil plane 
phasors. The Gaussian pdf of the image plane phasors is used 
to compute the pdf of the intensity, p{I). 

The form of p(I) follows from steps involving a change 
of variables from the phasor components. Originally derived 
by Goodman ( 1975), the resulting function is the modified 
Rician (MR) distribution, 

1 / I+IA flJTTA 
PMR(/)=-exp(^-— j/o^^j, (8) 

where /o(x) denotes the zero-order modified Bessel function 
of the first kind. For 4 = 0, the distribution becomes the ex- 
ponential statistics of 'pure' speckle, of primary interest in 
speckle interferometry. 

Intensity fluctuations play a critical rol e in the signal-to- 
noise ratio of direct imaging observations. IGoodmanl Jl975h 
showed that the mean and variance of / are 

^lI=k.+h (9) 

aj = lj + 2IJs. (10) 
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From equation (I1Q> . we see that the variance arising from pure 
speckle (the first term on the right-hand side of the expression) 
is augmented by the presence of the coherent wave. 

We have expressions describing the statistics of image 
plane intensity in an idealized absence of photon and de- 
tector noise. The effect of photon nois e on equations (Et 
and (I1Q> has been studied elsewhere JCana les & Ca gigall 
fr999aic .Aime & Soummed I200I ISoummer & Aimd 12004 
The additi on of in coherent light (e.g. from a planet) is ad- 
dressed by AS04. Detector read noise can be a significant 
factor in short-exposure images. We can model read noise as 
the sum of a Gaussian random variable with the actual speckle 
intensity. We take the read noise distribution to be zero-mean 
with standard deviation (Tm. Denoting this the MRG model, 
the pdf for the resulting observable is broader than the MR 
case, and is given by the convolution of equation with a 
Gaussian pdf, 

Pmrc(I) = Pmr(I)*Pc(I)- (11) 

3. OBSERVATIONS 
3.1. Measurements and Data Processing 

On 2004 July 31, we observed the star Sadalmelik (a Aqr, 
G2Ib), chosen for its apparent brightness in the K band (niK = 
0.59, 2.5 X 10^ photons m"^ fim~^ s~'). We obtained compen- 
sated short-exposure images under good weather conditions 
using the AO system of UCO/Lick Obser vatory's 3-mete r 
Shane Telescope on Mt. Hamilton ( Bauman et alJ I2OO2I) . 
There are 61 actuators on the deformable mirror used for 
wavefront control. The AO system was run at 500 Hz in a 
closed feedback loop. 

Rapid sequences of short-exposure images were obtained 
with IRCAL, the infrared camera which images the AO- 
corrected light onto a 256-bv -256 Rockwell Semiconductor 
PICNIC array (iLlovd etal .' 2000). Detector readout is driven 
by a digital signal processor within an Astronomical Read- 
out Cameras, Inc. electronics system (Leach & Low 200(3). 
The observations were made through a narrow-band Bij fil- 
ter (2.16 /im, AA/A ~ 1%), where the detector is Nyquist 
sampled by 75.6 mas pixels. 

We developed an array readout mode optimized for mini- 
mal delay between exposures in order to accurately measure 
the spatial and temporal distribution of stellar speckle pat- 
terns. Typical mode parameters were for correlated-double- 
sampling (CDS) a 32x32 pixel (2"4 x 2."4) subaiTay with 
5 ms integration times. With these settings, exposures were 
obtained every Af = 14.5 ms, and in 13 seconds a sequence 
of 900 images could be obtained before writing to disk. As 
demonstrated below, time spent in overhead between expo- 
sures (pixel 'dead time', about 9.5 ms) is unimportant, as it is 
shorter than the speckle decorrelation time (fdead < Tc)- The 
rms CDS read noise ai-„ is measured to be 45 e~ in this mode. 
Photon noise fluctuations do not strongly affect the statistics 
of these data. 

The image sequences were processed before analysis, in- 
cluding bias-subtraction and flat-fielding. Tip/tilt error affects 
image plane intensity statistics, a fact illustrated by consid- 
ering the extended Marechal approximation: the Strehl ratio 
(and by extension, total coherent intensity) decreases expo- 
nentially with the phase variance. We removed the residual 
tip/tilt error, reducing phase variance, to maximize the range 
of Ic- After a centroid measurement, each image was shifted 
and resampled with a damped-sinc interpolator The observed 
2-D rms tip/tilt error was .43 pixels (32 mas). 



In section |4] we will examine the distribution function for 
speckle intensities at each pixel using the measured time se- 
ries. Examples of these time series for three locations in the 
image plane are shown in Figure[2 The locations were chosen 
to highlight the diversity in intensity distribution as a func- 
tion of location in the image plane. Large deviations from 
the mean intensity are visible at location (a), while the data 
in location (c), close to the PSF core, do not show such large 
excursions. 

3.2. Temporal Characteristics 

We wish to analyze our measurements in light of the prob- 
ability distribution functions for the intensi ty ra ndom vari- 
able. The statistical tests we will employ in ^4. II require that 
the tested data be drawn from independent random variables. 
However, we observe the intensity as a realization of a tempo- 
ral random process, and correlations exist between intensities 
at different times. In order to perform tests concerning the 
distribution of independent values of intensity, we must de- 
termine the timescale over which samples of I(t) can be con- 
sidered independent. We assume temporal stationarity for the 
duration of the observations, so that we may take the statisti- 
cal properties to be constant in time. 

At a single spatial location in the tip/tilt-removed image 
sequence, the number of independent samples can be char- 
acterized by the speckle decorrelation time, rdx). Previ- 
ous definitions of Tc have been motivated by the the desire 
to compute optimal exposure times in speckle interferom- 
etry ( Roddier. Gilli. & Lund i_9_82i and references therein). 
The decorrelation time is derived from the temporal auto- 
covariance function C/(t). It can be the e-folding va lue of 
this fu nctio n, as in Scaddan & Walker ( 1978), Roddi er et alJ 
(1982'), and'Ver ninet alJ (I199L) . It can also be defined by the 
equivalent width, 

/oo 
Cj(T)dT, (12) 
00 

following lAime et all J1986I) . who found their empirical 
C/(r) was fit by the sum of two Lorentzian functions. 
The equivalent-width definition assumes no prior knowledge 
about the shape of the autocovariance function, and has 
the advantage of using information contained in the func- 
tion's tail. Decorrelation times Tc computed from equa- 
tion J12t are generally longer than those computed by e- 
folding (Aime et al. 1986). This definition also has the prop- 
erty that the variance of the intensity integrated over time 
T ^ Tc is given by ajT /re- 
in practice, we compute r^. for each pixel, approximating 
the integral as the sum of an estimator for the autocovariance 
sequence {acvs) over the interval ±35 samples (±0.5 s). We 
restrict the summation domain because the acvs estimate be- 
comes more uncertain with increasing lag t. Estimating the 
error in is complicated by the large correlations between 
terms in the acvs estimate, and it suffices for the analysis in 
^4. II to neglect its calculation. Limiting the summation do- 
main in the calculation of r^. generally results in an underes- 
timate of the true decorrelation time. Finally, note that the 
effect of read noise is not removed from the acvs in the 
calculation. In this case, Tc will be lower than the speckle- 
only decorrelation timescale in regions of significant detector 
noise. 

Figure|2]shows the acvs estimates for the selected locations 
from Figure^] along with the given by the aforementioned 
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Fig . 1 . — The long-exposure image (left, log scale), and a plot of observed intensity as a function of time for three selected pixels (right), the encircled locations 
(a), (h), and (c). The locations were chosen to highlight the diversity of speckle intensity distribution as a function of image plane position, moving from the PSF 
halo to the first Airy ring. For each pixel, 900 measurements are shown at a sampling period of 14.5 ms, which consists of 5 ms integration and 9.5 ms of dead 
time. To illustrate differences in speckle fluctuations with field position, the intensity for each location is normalized by its sample mean, which was 54, 200, and 
1800 respectively for (a), (b), and (c). The scale of fluctuations relative to the mean intensity decreases from locations (a) to (c). Read noise and decoirelation 
time are shown as scale bars. The read noise at location (c) is neghgible. 

procedure. As shown in the inset of Figure |2j the are ob- 
served to vary with field position x. Analysis of the spatial de- 
pendence of speckles' temporal decorrelation is deferred for 
future work. 

The measured decorrelation times are generally many times 
the sampling rate At, which ensures that we may disregard 
the effect of finite integration time and pixel dead time on the 
analysis of the observed intensity distributions. The median 
decorrelation time was 175 ms, and the first and third quar- 
tiles of the measured r^Cx) distribution are 95 ms and 270 ms, 
respectively . 

As in Aime et al. ( 1986), the forms of the acvs estimates 
in Figure |2l are qualitatively consistent with the sum of two 
Lorentzian functions. A unique decorrelation timescale char- 
acterizes the width of each component Lorentzian. For lo- 
cations (a)-(c), the fast timescale varies between 20 and 80 
milliseconds while the slow timescale lies between .8 and 1.2 
seconds. 

4. ANALYSIS 
4. 1 . Testing the Distribution 

The observed distribution of speckle intensities can act as 
a consistency check against the modified Rician distribution, 
Pmr{I) defined in equation (|8}. Neglecting photon count- 
ing fluctuations, we expect the observed intensities to be the 
sum of a modified Rician random variable and an independent 
Gaussian random variable for the read noise, with distribution 
Pmrg(I) given by equation il 1> . We adopt as a null hypothesis 
that the distribution of intensities follows the model distribu- 
tion, with parameters and I, estimated from the data (where 
the caret denotes an estimate). 

For each pixel location, the observed sequence of images 
gives a measure of the intensity as a function of time. Data 
from five consecutive 900-exposure sequences were com- 



FlG. 2. — Autocovariance sequence (acvs) estimates for the time series 
intensity data from pixels (a)-(c) in Figure Q The abscissae of the filled 
circles denote the decorrelation times Tc computed from the acvs estimates. 
The Tc characterize the number of independent speckle realizations in a given 
time period, and are observed to vary systematically with field position as 
illustrated by the inset map. Sampling effects are minimized with the short 
sampling period of 14.5 ms. 



bined into the same histogram. Section l372l describes tempo- 
ral correlations in the observed sequences, which arise from 
correlations in the perturbations a(t,u). To mitigate the ef- 
fect of temporally correlated data in hypothesis tests which 
assume independent measurements, we formed the histogram 
after decimating each intensity sequence at intervals of the 
pixel's decorrelation time. At' = Tc- The data in between these 
intervals were discarded, as were values falling below zero. 
Examples of the resulting histograms are shown in Figure |3 
using the data for locations (a)-(c) from Figures^and^ 

An iterative least-squares fit of the model distribution to the 
histogram estimates the parameters and Is, holding the read 
noise (Tm fixed. We perform this fit at each location in the im- 
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Fig. 3. — Histograms of intensity with best-fit model distribution function 
(dashed curve) for the data of locations (fl)-(c) shown in Figures □andlH Also 
shown are l-cr en'ors for the model distribution. As in Figure II I intensity 
is normalized by the sample mean. Inset are the number of samples in this 
histogram (after resampling the original timeseries by A?' = Tc), and the ratio 
of parameters Ic and /, characterizing the model distribution. 

age. Best-fit models for distributions of locations (fl)-(c) are 
also shown in Figure|3] We choose the least-squares approach 
to parameter estimation for simplicity, though one might de- 
velop maximum-likelihood estimators. We briefly discuss po- 
tential bias in our estimators below. 

The residuals of the fit are used in a test to determine 
if the hypothetical model distribution function should be re- 
jected. The number of degrees of freedom (nominally the 
number of histogram bins) is reduced to account for the es- 
timation of distribution parameters from the data. We can 
calculate the p-value from the distribution and reject the 
hypothesized distribution if this measure falls below the sig- 
nificance level a (taken to be 5%). 

In contrast to a generic test, we can also use a hypothesis 
test based on the empirical cumulative distribution function 
like the Kolmogorov-Smirnov (K-S) test. While both tests 
determine a p-value for rejecting hypotheses, we expect the 
K-S test to be make fewer Type II errors (failing to reject 
false null hypotheses) because it avoids binning data. Dis- 
advantages to the K-S test include a relative insensitivity to 
distribution tails and a much larger computation requirement; 
as we are using estimated parameters in the hypothesized dis- 
tribution, it is necessary to determine the K-S test statistic's 
distribution function to find the p-value. This can be done 
with a bootstrap Monte Carlo method, where the parameters 
Ic and Is (estimated from the observations) are used to simu- 
late a random variable following the model distribution. For 
each trial / in the Monte Carlo simulation, distribution param- 
eters are estimated from the simulated data (via the preceding 
least-squares fit procedure) and the K-S test statistic D, is cal- 
culated. After the completion of 1,000 trials, the distribution 
of Di is used to estimate the /9-value. We restrict the use of this 
form of the K-S test to example pixel locations (a)-(c) because 
the requisite simulations cannot be calculated in a timely fash- 
ion given our computational resources. However, the test 
is applicable to the entire field. 

The results of the hypothesis tests for the example locations 
are shown in Table [2 We see that neither the nor the K-S 
tests reject the hypothetical MRG model for these data. The 
table presents similar test results for the MR model, as well 
as a Gaussian (G) with free parameters /i and a. We did not 
test a separate Poisson model because the data are not in the 
photon counting regime. As read noise significantly affects 



TABLE 1 
Hypothesis Test Results 



Location 


Model 


X"/«dof 


test 


p- value 

K-S test 


(a) 


G 


1.61 


0.051 


0.26 ± 


0.03 


(a) 


MR 


2.10 


0.0087 


0.001 ± 


0.001 


(a) 


MRG 


0.90 


0.46 


0.11 ± 


0.02 


(b) 


G 


2.51 


0.0011 


0.015 ± 


0.007 


(b) 


MR 


0.51 


0.87 


0.001 ± 


0.001 


(b) 


MRG 


0.55 


0.83 


0.10 ± 


0.02 


(c) 


G 


1.44 


0.096 


0.17 ± 


0.02 


(c) 


MR 


0.85 


0.51 


0.52 ± 


0.03 


(c) 


MRG 


0.85 


0.51 


0.54 ± 


0.03 



Note. — Results of x' and Kolmogorov-Smirnov hypothesis tests for fits 
to intensity data. G, MR, and MRG refer to the Gaussian, modified Rician, 
and modified Rician plus a (zero-mean, a = (Tm) Gaussian models for inten- 
sity distribution. Instances of distribution rejection are indicated in italics, 
and occur when the p-value falls below the significance level a = .05. The 
mean and variance of the G model and parameters Ic, h of the MR and MRG 
models are estimated from the data, and jijof is the number of histogram bins 
minus the number of estimated parameters. The quoted uncertainty in the 
K-S test p- values coiTesponds to 95% confidence. 

the observed distribution of location (a), the model distribu- 
tions consistent with Gaussian read noise provide the best fits. 
However, the pure Gaussian model G is only marginally sig- 
nificant for locations (a) and (c), while the asymmetric dis- 
tribution at {b) ensures rejection. The MRG model is not re- 
jected by either test at all three locations. 

When examining the results of the x^ hypothesis test over 
the image, we find that the MRG distribution hypothesis is 
not rejected for 90% of the pixels. This is slightly lower 
than the expected rate of (1 - a) = 95%, which is indicative 
of marginally significant processes which alter the intensity 
distribution. 

Several processes may bias the results of these tests. Read 
noise acts to broaden the density function, which limits the 
accuracy and precision of the Ic and estimates in regions 
of low signal. Additionally, we have chosen to neglect photon 
noise, as generally ^ (/). Finite pixel size also plays a role 
in altering the distribution from the ideal case, as the pixels' 
area is a significant fraction of a speckle. The spa tially inte- 
grated speckle distribution has also been studied bv lGoodmanI 
( 1975, 1985), though we do not consider that case here. An- 
other caveat arises from the non-stationarity of the process 
governing the wavefront distortion term fl(u,0 in the com- 
pensated imaging sys tem, as this te rm is driven by the at- 
mospheric turbulence llRoddieill98lL and references therein). 
These fluctuations can bias the estimates Ic and Ic and there- 
fore reduce goodness-of-fit. Finally, we expect second-order 
correlations in the sampled intensities arising from under- 
estimates to cause Type I errors, an increase in the test rejec- 
tion rate from a when the hypothesis is true. 

4.2. Coherent PSF Extraction 
The parameter estimation procedure above provides the 
spatial distributions of (x) and /.s(x), the sum of which forms 
the long-exposure image via equation 0. Figurel^shows the 
/f(x) and /5(x) extracted from the observations. As we obtain 
uncertainties of the estimates for these quantities via least- 
squares fitting, it is possible to make a crude calculation for 
the S/N in the estimated parameters. Neglecting systematic 
biases, the median S/N over all pixels are 7.6 and 4.8 for /^(x) 



6 




<I> I,(x) I,(x) 

Fig. 4. — Decomposition of the long-exposure PSF ~ (/(x)), reprised from 
Figurefn into /c(x) and /,(x) as in equation j4). By definition, tlie quantity 
/c(x) is proportional to the static PSF. Images are displayed with the same 
logarithmic scale. 

and /i(x), respectively, in data gathered in 1 minute of observ- 
ing time. It may be possible, through maximum-likelihood 
analysis, to develop improved estimators which avoid binning 
data before parameter estimation. 

While /f(x) is a difficult quantity to measure, as it is pro- 
portional to the PSF of a system free of atmospheric turbu- 
lence, its knowledge can be valuable. It may have utility 
in calibrating methods which reconstr uct the long-exposur e 
PSF (T > Tc). With the technique of lVeran et all jl997aft . 
data from the AO system are used to calculate a blurring ker- 
nel. This kernel is convolved with a static PSF (cx Ic) to es- 
timate the on-axis long-exposure PSF. This static PSF con- 
tains non-common path aberrations, and in standard practice 
is calibrated by an artificial point s ource or observations of 
a reference star SV€mn et alJll997bl) . A drawback in the ar- 
tificial source case is that aberrations with high spatial fre- 
quency arising from the primary mirror are not sensed. When 
using a reference star, calibrating the static PSF requires an 
accurate construction of the blurring kernel during the refer- 
ence observation. Obtaining /^(x) by estimating the parame- 
ters of the speckle intensity distribution function (as in ^4.1> 
can be used to calibrate the PSF reconstruction in a manner 
independent of a blurring kernel estimate, though errors such 
as those arising from temporal non-stationarity may limit the 
accuracy. Finally, we note that the reference star calibration 
method presented in Veran et al. ( 1997b) and the speckle dis- 
tribution function method above require that the calibration 
target be a point source. If incoherent light is added as from 
a companion or extended emission, the speckle distribution 
will be altered from the MR. An analytical expression for the 
intensity distribution function in this case is given bv .AS04. 

5. CONCLUSIONS 

We have shown observations of speckle intensity over 
timescales of ^1 min are consistent with modified Rician 
statistics through and Kolmogorov-Smirnov goodness-of- 
fit tests. In high-contrast imaging for planet detection, the 



variance of the integrated intensity at the location of the planet 
sets the noise level against which the signal must be detected. 
Central to understanding this noise level are the variance of 
intensity <jj and the speckle decorrelation time t^. AS04 high- 
light the efficacy of coronagraphs in reducing crj by reducing 
terms containing Ic- For example, suppose Ic = Is= 1- From 
equation ( II Ot . aj = 3. With a coronagraph, Ic = and aj = 1, 
a reduction by a factor of 3 in variance. However, a complete 
understanding of speckle fluctuations in coronagraphic imag- 
ing experiments will require characterization of t^, as well as 
understanding the effects of non-stationarity in Ic and /, due 
to changes in the atmosphere and quasi-static wavefront errors 
in the optical system. Analyses of observed intensity distribu- 
tions over timescales longer than 1 min are left to future work. 

We note that the dark speckle imaging technique relies on 
analysis of t he obs erved intensity distribution at each location 
X (Labevrid ri995t iBoccaletti et al. 1998), and has tradition- 
ally assumed exponential statistics for intensity. Here we have 
demonstrated that probability distributions related to the mod- 
ified Rician are applicable in the case of moderate wavefront 
correction and thus may be incorporated into this technique. 
AS04 give speckle distributions for additive planet light and 
photon-counting conditions which can be used to this end. 

Additionally, we have presented a technique where /^(x) 
and /5(x) are estimated from the observed distributions of in- 
tensity. Knowledge of the unaberrated PSF, proportional to 
/f (x), can b e applied to the cali bration of PSF reconstruction 
algorithms ( I Veran et alJll997al) . This application may spur 
the further development of estimators for Ic and Is given an 
observed speckle intensity distribution. 
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